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Summary 

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) is modified 
to allow the calculation of turbulent flows. This is accomplished using the Cebeci- Smith 
and Baldwin-Lomax eddy-viscosity models in conjunction with the thin-layer Navier-Stokes 
options of the program. Turbulent calculations can be performed for both perfect-gas and 
equilibrium flows. However, a requirement of the models is that the flow be attached. It is 
seen that for slender bodies, adequate resolution of the boundary-layer gradients may require 
more cells in the normal direction than a laminar solution, even when grid stretching is 
employed. Results for axisymmetric and three-dimensional flows are presented. Comparison 
with experimental data and other numerical results reveal generally good agreement, except 
in the regions of detached flow. 


Introduction 

With the continuing interest in hypersonics, and the associated desire to model more 
complex phenomena of the flowfields about bodies in hypersonic flight, numerical solvers 
of the full Navier-Stokes equations are being employed with increasing regularity. Ongoing 
advancements in numerical algorithms and computer hardware have made this feasible. Spe- 
cific applications include full-flowfield calculations about complex configurations such as the 
Space Shuttle 1)2 and National Aero-Space Plane (NASP) 3 ’ 4 , with their associated high 
altitude, high speed environments. Development of the various numerical models necessary 
to simulate the physical processes in these environments is an important part of this capa- 
bility. At high Reynolds number conditions, turbulent flow is one process which becomes 
important, and its numerical simulation is the subject of this work. 

Solving the full Navier-Stokes equations is required to accurately predict most compli- 
cated continuum flowfields. Typically the unsteady form of the Navier-Stokes equations are 
solved using a time-marching procedure. This allows any elliptic behavior to be properly 
modeled, but at a price: Navier-Stokes solutions are significantly more computation-intensive 
than solutions from less detailed methods which do not involve time relaxation procedures. 
The thin-layer approximation to the Navier-Stokes equations may be employed (to eliminate 


the viscous terms containing derivatives in the streamwise direction), but the solution is 
still quite expensive. Still, for complex flowfields with separation and large subsonic regions, 
there is little choice but to solve the full or thin-layer Navier-Stokes equations. 

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) is a three- 
dimensional finite- volume Navier-Stokes solver developed by Gnoffo 5-8 . LAURA was ini- 
tially applied to blunt-body flows, but has been used with slender bodies 9> 10 as well. This 
paper serves as a follow-up to reference 10 w T hich used LAURA for aerothermodynamic pre- 
dictions over slender vehicles assuming laminar, perfect-gas flow. In this work, turbulent 
flows over slender bodies will be modeled. 

In various efforts 11-13 eddy-viscosity models have been mated with existing flowfield 
solvers. The Cebeci-Smith ^ and Baldwin-Lomax ^ eddy-viscosity models are employed 
here. Such algebraic (or zero-equation) models are less complicated than more exact ap- 
proaches, such as the Johnson-King 18 and two-equation 13 models, and as a result are more 
computationally efficient (although theoretically less accurate). The thin-layer Navier-Stokes 
option of LAURA is exercised in this study. For implementation in LAURA, several issues 
must be addressed: 1) the actual modeling of the turbulence, 2) the resolution of near-body 
gradients, and 3) the transition from laminar to turbulent flow. 

In this study, both perfect-gas and equilibrium flow conditions are considered. Calcu- 
lations are performed for sphere-cones and a generic transatmospheric vehicle (TAV). The 
LAURA heat-transfer predictions are compared with experimental data, as well as other nu- 
merical calculations, in order to evaluate the present results. In the sections that follow, brief 
discussions of the LAURA computational method, along with the modifications employed 
for turbulent calculations, are given. A discussion of the results and some conclusions from 
this study follow. 


List of Symbols 


a 

Speed of sound 

C v 

Specific heat at constant pressure 


Unit vector in s-direction 


Unit vector in ^-direction 

h 

Static enthalpy 

H 

Total enthalpy; H = h -f V 2 / 2 

k 

Thermal conductivity 

L 

Reference length of body 

M 

Mach number 

n 

Normal distance from the body 

A n 

Height or thickness of cell (in n-direction) 

n + 

Normal coordinate parameter for turbulent flow 

Pr 

Prandtl number; Pr = f. iC p /k 

q 

Heat transfer rate 

Re ce n 

Cell Reynolds number; Re ce u — paAn/fi 

s 

Surface distance along the body, measured from its nose 

T 

Temperature 

U, V, w 

Velocity components in the s-, n-, ^-directions, respectively 
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V 

Vr, V„ V Mt 

x, y , z 
6 
6 * 
e + 

V 
P 


u 


Total velocity; V 2 = V 2 + V 2 + V 2 = u 2 + v 2 + w 2 

Velocity components in the x-, y-, ^-directions, respectively 

Cartesian coordinates with origin at body nose 

Value of 7i at boundary-layer edge 

Boundary-layer displacement thickness 

Ratio of eddy viscosity to laminar viscosity; e + = n t /p 

Dynamic viscosity 

Density 

Angular coordinate in the circumferential direction around the body 

Shear stress 

Vorticity 


Subscripts 

e Boundary-layer edge 

t Turbulent value 

W Wall value 

oo Freestream condition 


Computational Method 

LAURA is a finite-volume formulation of the integral form of the Navier-Stokes equa- 
tions. A second-order accurate, symmetric total-variation-diminishing (STVD) scheme 19 is 
used in conjunction with the upwinded differencing of the discretized equations. At each 
cell face, Roe’s averaging 20 defines the flowfield variables based on values from the adjacent 
cells. The unsteady governing equations are driven to a steady-state solution through a 
time relaxation procedure involving global sweeps through the computational domain. In 
the relaxation scheme, dependent variables at a cell center are treated implicitly, using the 
latest available information from adjacent cells. Thus, for perfect-gas flow, the relaxation 
simply requires the inversion of a 5 x 5 matrix at each cell center. 

During the relaxation process, the grid is periodically adapted in the body-normal di- 
rection so that the grid can be tailored to the emerging solution. This process of customizing 
the mesh involves clustering the cells near the body surface for accurate heat-transfer cal- 
culations. In addition, the outer boundary of a converged grid is parallel to the captured 
bow-shock wave. A later section provides more details of this process. 

Typically, calculations begin with a grid which has only one-fourth the total number of 
cells in the normal direction. The use of this coarse grid allows even a poor initial guess for 
the mesh to quickly align itself with the developing bow shock with repeated adaptions of 
the normal distribution. The density of cells in the normal direction is later increased by 
factors of two during the relaxation sweeps. This process promotes convergence and stability 
of the solution. 

For future reference, the STVD limiter given by equation (3.8b) of reference 7 was used 
for the results presented later. Further, for the present calculations, an eigenvalue limiter of 
0.3 was used in all directions. For the normal direction, the limiter is scaled according to 
the cell aspect ratio so that the resultant limiter approaches zero near the body. As a result, 
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the viscous terms containing derivatives in the streamwise direction), but the solution is 
still quite expensive. Still, for complex flowfields with separation and large subsonic regions, 
there is little choice but to solve the full or thin-layer Navier-Stokes equations. 

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) is a three- 
dimensional finite- volume Navier-Stokes solver developed by Gnoffo 5-8 . LAURA was ini- 
tially applied to blunt-body flows, but has been used with slender bodies 9> 10 as well. This 
paper serves as a follow-up to reference 10 which used LAURA for aerothermodynamic pre- 
dictions over slender vehicles assuming laminar, perfect-gas flow. In this work, turbulent 
flows over slender bodies will be modeled. 

In various efforts 11-13 eddy-viscosity models have been mated with existing flowfield 
solvers. The Cebeci-Smith 14-16 and Baldwin-Lomax 17 eddy-viscosity models are employed 
here. Such algebraic (or zero-equation) models are less complicated than more exact ap- 
proaches, such as the Johnson-King 18 and two-equation 13 models, and as a result are more 
computationally efficient (although theoretically less accurate). The thin-layer Navier-Stokes 
option of LAURA is exercised in this study. For implementation in LAURA, several issues 
must be addressed: 1) the actual modeling of the turbulence, 2) the resolution of near-body 
gradients, and 3) the transition from laminar to turbulent flow. 

In this study, both perfect-gas and equilibrium flow conditions are considered. Calcu- 
lations are performed for sphere-cones and a generic transatmospheric vehicle (TAV). The 
LAURA heat-transfer predictions are compared with experimental data, as well as other nu- 
merical calculations, in order to evaluate the present results. In the sections that follow, brief 
discussions of the LAURA computational method, along with the modifications employed 
for turbulent calculations, are given. A discussion of the results and some conclusions from 
this study follow. 


List of Symbols 


a 

Speed of sound 

c p 

Specific heat at constant pressure 


Unit vector in s-direction 

e<t> 

Unit vector in (^-direction 

h 

Static enthalpy 

H 

Total enthalpy; H = h + V 2 /2 

k 

Thermal conductivity 

L 

Reference length of body 

M 

Mach number 

n 

Normal distance from the body 

An 

Height or thickness of cell (in n-direction) 

n + 

Normal coordinate parameter for turbulent flow 

Pr 

Prandtl number; Pr= pC p /k 

? 

Heat transfer rate 

Pt-cell 

Cell Reynolds number; Re ce ii = paAn/fi 

S 

Surface distance along the body, measured from its nose 

T 

Temperature 

U , V, w 

Velocity components in the s-, n-, ^-directions, respectively 



Conductivity is also affected by turbulence. For laminar flow, 

„ Jl?!l 

q ~~C ? dn 


Accounting for the turbulent contribution gives 

(k + k t ) dh 
q= C p 


( 8 ) 

( 9 ) 


where kt is the eddy thermal conductivity. By definition, 


Pr = 


CjP_ 

k 


so that 


k /x 

~C~ V ~ ~Pr 


and 


k t lit 

C p ~ Prt 


Thus, 

A + i± = AiJfL = iL(i + ftifL] 

C p C p Pr Pr t Pr\ /x Pr t ) 

While in general the value of Prt varies throughout the shock layer, a constant value of 0.9 
is commonly used. 

Defining the ratio of eddy viscosity to laminar viscosity as 


£ 


+ 


P± 


( 10 ) 


gives the final form of the modifications to the governing equations: 






~( 
Pr V 


1 + e H 


Pr\ 

Pn) 


Note that for laminar flow, fi t = 0 and k t = 0. 

The Cebeci-Smith and Baldwin-Lomax algebraic models are two-layer eddy-viscosity 
formulations which are applicable to attached flows. The inner-layer value, ef , is based on 
Prandtl’s mixing length concept 22 . For Cebeci-Smith, the outer layer value, e* , is given 
by the Clauser-Klebanoff 23, 24 expression. The Baldwin-Lomax model uses an equivalent 
expression for £^, but differs in its determination of the length scale. For both models, the 
inner-layer expression is used from the wall outward until ef > ef , thus forming a composite 
eddy viscosity. Details of each of these models are given below. 


Cebeci-Smith 

As mentioned above, the inner-layer eddy viscosity is based on Prandtl’s mixing length 


theory: 



dU 

dn 


( 11 ) 
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The mixing length for the inner layer according to Van Driest 22 is 

l = K v n J'l — exp [—n + /A + y^ 


(12) 


where the von Karman constant is 


and the normal coordinate parameter is 


K v = 0.4 


+ n P 
n = — 


1 1/2 


(13) 


In reference 25, the damping factor is defined as 


A + 


-»© 


where the local shear stress is 


and 


r = p (l + 


- 1/2 


dU 


dn 


I'll) — fJ'h 


dU 


dn 


(14) 


(15) 


(16) 


The outer layer eddy viscosity is approximated by the Clauser-Klebanoff 23, 24 expression: 


et = 0.0168 U r 6 * 




The displacement thickness is 




(17) 


(18) 


w’here 6 is the value of n at the boundary-layer edge. By definition, at the boundary-layer 
edge, 

H 

uZ 1 

Numerically, this can be approximated as the grid point where 

H 


H a 


> 0.995 


The normal intermittency factor is 


7i,n 


1 + 5.5 


(i) 


“1 


(19) 


6 


( 12 ) 


The mixing length for the inner layer according to Van Driest 22 is 

l = K v n 1 — exp (— n + /.A + )j 
where the von Karman constant is 


K v = 0.4 


and the normal coordinate parameter is 


+ n P 

n = — 


1 1/2 


P L P 

In reference 25, the damping factor is defined as 


(13) 


A + 


-»(ia 


where the local shear stress is 


and 


T = p (l + £ + ) 


- 1/2 


dU 


dn 


T w Ph 


dU 


dn 


(14) 


(15) 


(16) 


The outer layer eddy viscosity is approximated by the Clauser-Klebanoff 23 ’ 24 expression: 


£+ = 0.0168 Ue 6* 


P 7«> 

P 


The displacement thickness is 




dn 


(U) 


(18) 


where 6 is the value of n at the boundary-layer edge. By definition, at the boundary-layer 
edge, 


H 

HZ 


1 


Numerically, this can be approximated as the grid point where 

H 


H n 


> 0.995 


The normal intermittency factor is 


7«,n 


1+5 


■>© 


-I “I 


(19) 




Baldwin-Lomax 


The expression for the inner layer eddy viscosity in the Baldwin-Lomax model is similar 
to that used in the Cebeci-Smith model: 



( 20 ) 


The mixing length, l, is again given by equation 12 in conjunction with equations 13 
through 16. The magnitude of the vorticity, |w|, is found from 


M = 


dV 2 dV v 


12 


dy 


dz 


+ 


dV z dV x 
dx dz 


+ 


dVy dV z 

dx dy 


1 2 'I 2 


dw dv 
[ dn d(f> J 


+ 


dw du 
ds d(j) 


1 2 


+ 


dv du 
ds dn 


I 2 | 2 


For thin-layer Navier-Stokes, this becomes 


u> = 


[ dw 
dn 


1 2 


+ 


du 

dn 


12^2 


( 21 ) 


( 22 ) 


The outer layer eddy viscosity is approximated by the Clauser-Klebanoff : 

£+ = 0.0168 C C p F, - P 7,> 


max Tl-max 




expression: 

(23) 


where n — n max is the location of the maximum value, F max , of the vorticity function: 


F = 


n \u\ 


exp n + /A + )J 


(24) 


Although a Mach number dependency has been suggested for Cep (see ref. 25), in the 
results presented later, a constant value of Ccp = 1.6 as given in reference 17 yields the best 
agreement. Klebanoff’s intermittency factor is 


7«,n — 


1+5.5 


Ckleb n \ 6 

W max ' 


(25) 


where Ckleb = 0.3. 

The boundary-layer thickness does not appear in the Baldwin-Lomax model. The length 
scale of the outer layer is instead based on F max , the maximum of the vorticity function. Since 
the solution is discretized, the position and value of this maximum can be “smeared”. In 
an effort to minimize this behavior, the discrete values in the vicinity of the maximum are 
curve-fit. using an “overlapping-parabola” technique. The result of this process can be seen 
in the sample streamwdse distributions shown in the figures which follow. Figure 1 shows 
that the process has little effect on F max itself. However, figure 2 shows that the position of 
this maximum is a more sensitive quantity. Note the “jump” at s fa 4 which is present when 
the parabolic blending is not used. Since it is the product of these terms, F max n max , which 
appears in equation 23, the calculated eddy viscosity will reflect such “jumps”. Figure 3 
illustrates this effect on the near-body ^-distributions in the vicinity of a jump in n max . 
Note that the streamwise variation in the £ + -profiles is much smoother when parabolic- 
blending is employed. 


7 



where Re ce u }W is the desired cell Reynolds number at the body and n k is the distance from 
the body to the outer grid boundary. The value of uk is “lagged” from the previous grid. In 
past investigations 10 , a wall value of Re ce u iW = 0(1) at the wall has proved to be adequate 
for laminar heating rates. 

The objective is to adequately resolve the boundary-layer gradients without overly- 
resolving the outer shock layer where gradients are small. Thus, grid stretching is used near 
the body, in conjunction with even spacing in the outer region of the grid. The spacing from 
these two regions are constrained to match at their interface. The fraction of the total cells 
used in the stretched region is given by 

/ l 32 \ 

F s tr = max ^ 2 ’ 1 ~ A 7 / ( 31 ) 


where K is the total number of cells in the normal direction. The number of stretched cells 
for this normal line is 

Kstr — F s t r K (32) 

The recursion formula for the nondimensional height of these K str cells, using a sinusoidal 
distribution, is 

(k — 1) 7T j 


A h k = < 1 + f stT sin 


A sir — 1 


An*_i 


where 


fstr — 


F, 


sir 


v sir 


i 


(33) 


(34) 


lAni J 

The total nondimensional distance from the body to the outer edge of any cell k is found 
from 

k 

= (35) 

/=i 

This distribution is then normalized to yield 0 < h < 1. 

In this investigation, it was found that the above normal distribution did not always 
adequately resolve the turbulent boundary layer. Apparently, the cell Reynolds number 
at the wall must be more tightly matched to the target value of Re ce n >w = 1. To achieve 
this, the value of f 3tr is determined iteratively for each stack of cells, with the value from 
equation 34 serving as an initial guess. The resultant iterated value of / 5<r yields fix ~ 1 so 
that Re C€ u,w ^ 1 at the wall. 

Another quantity useful in gauging adequate grid resolution for turbulent flows is the 
wall value of the normal coordinate parameter (n+). As a rule of thumb, its value should 
not be larger than (9(0.1) for proper resolution of the turbulent boundary-layer gradients. 
Using a target value of = 0.1, this yields another expression for the thickness of the cell 
adjacent to the body: 

(36) 


nt n\ 


A hr = 

nj no- 


where, from equation 13, 


+ _ Til pi 


V 1 


1 1/2 


nv = 


Pi. 


9 


PRESIDING PAGE. BLANK NOT FILMED 


Note that n, is the distance from the body to the center of cell 1, so that n, = 2 n,. 

Both equations 30 and 36 provide values for An,, and so the smallest of these two values 
is used to define the thickness of any cell adjacent to the body. With the proper stretching 
determined, dimensionality is returned to the distribution by 

n k = y h k (37) 

where n s is the distance from the body to the bow shock for this stack of cells (lagged from 
the previous grid). The fraction of the distance from the body to the outer grid boundary 
which lies between the bow shock and body is given by F s . A value of F s = 0.8 is used here 
so that with a converged grid, the distance across the shock layer is eighty percent of the 
total distance from the body to the outer (freestream) boundary. 


Implementation 

The turbulence models used here require a fairly developed shock layer for implementa- 
tion. Therefore, for all the cases presented herein, the LAURA code is run for laminar flow 
until the residual drops several orders of magnitude. At this point (before the laminar solu- 
tion is fully converged), the turbulence is “turned on” and the iteration process is continued. 
Not surprisingly, the switch from laminar to turbulent flow results initially in a jump in the 
residual. After a period of adjustment to the new governing equations, the residual again 
begins to steadily decrease. The rate of this decrease is less than for laminar flow, since 
the higher viscosity for turbulent flow serves as damping for the solution. Figure 4 gives a 
typical residual history. In addition to the rise with the switch from laminar to turbulent 
flow, “jumps” in the residual occur when the grid is adapted. 

The LAURA package has several features which enable it to handle very large jobs 
and utilize multiple processors when running on CRAY supercomputers. For instance, the 
computational domain can be divided into several ’’blocks” with information communicated 
across their boundaries. The work within a given block (or for the entire domain) may be 
allocated to several processors. Although the results are not included here, the turbulence 
programming has been successfully tested in conjunction with both of these features. 

Solving the governing equations for turbulent flow increases the required computational 
time by approximately 1.5 percent per iteration per grid cell as compared with laminar 
calculations. However, since the turbulent terms increase the damping, more iterations are 
required to drive a turbulent solution to the same level of convergence as a laminar solution. 
Storage requirements for the LAURA code are also increased slightly since values for the 
eddy viscosity and normal distance from the body must be saved for each cell in the grid. In 
addition, the distance along the body surface to each wall-bounded cell is required to define 
the transition region. 


Results and Discussion 

In an effort to verify proper implementation of the turbulence models within LAURA, 
several test cases were computed. The results are given below in order of increasing complex- 
ity. The first case is an axisymmetric, perfect-gas flowfield calculation which is compared 
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with experimental data. The second is an axisymmetric, equilibrium flow case which is com- 
pared with results from another numerical approach. Finally, results for perfect-gas flow 
over a three-dimensional body are compared with experimental data. The results presented 
here are restricted to heating-rate comparisons. Specifically, streamwise heating rates (nor- 
malized by a reference rate) are plotted as a function of axial distance from the body nose 
(nondimensionalized by a reference length). For the three-dimensional case, circumferential 
heating rates for selected stations are included, plotted versus spanwise distance from the 
symmetry plane. 

Perfect Gas 

The first case is a Mach 5 perfect-gas flowfield calculation over an 8 -deg sphere-cone 
at 0-deg angle of attack. Comparisons are made with Jackson’s experimental data 29 . The 
geometry has a length of L = 25.4cm (l(hn) and a nose radius of R n ose — 6.35 cm (2.5m). 
The freestream conditions are p ^ = 0.48 kg/m 3 and = 111.1 I\. Jackson’s data include 
body temperature measurements, so a variable wall temperature is specified. The transition 
point is specified to be s = 2.54cm (1 in). Calculations using both the Cebeci-Smith and 
Baldwin-Lomax models are presented. The computations were performed on a workstation 
with 41 cells used in the axial direction. 

The baseline LAURA code can readily accommodate increases in the number of cells in 
the normal ( k ) direction. However, these increases are limited to grid doubling. For example, 
if 64 cells do not yield a satisfactory solution, then the grid may be doubled to 128 cells. 
For this investigation, the doubling routine has been modified so that any increase in the 
number of cells in the normal direction may be performed easily. To reduce computational 
costs, LAURA starts the solution for each new number of cells from the previous solution, 
rather than restarting from conditions based on freestream values. 

The heating results of a grid refinement study (using the Cebeci-Smith model only) are 
presented in figure 5, where q Te j = 340 kWjm 2 (30 BTU/s/ft 2 ). The condition of Re C eii,w = 1 
is enforced for all solutions shown in the figure. Recall that the total number of cells (7\ ) 
dictates the value of I\ str , the number of cells in the near-wall region (see equations 31 
and 32). Apparently, for this slender body, a value of K str > 48 is required to obtain a 
grid converged prediction of the turbulent heating. If a value of K str < 48 is used, a grid 
which extends out far enough to encompass the bow shock while maintaining Re ct \\ }W = 1 
is overly stretched. To further illustrate this point, table 1 gives the maximum stretching 
(1 + fstr) associated with each normal cell distribution represented in figure 5. Although 
this maximum stretching does not occur at the wall (see equation 33), it does lie within 
the boundary layer. It is seen that for this slender body, the K = 96 solution is essentially 
the same as the K = 128 solution. This is significant since the 128-cell solution requires 
approximately one-third more time to run than the 96-cell result. One can infer from these 
results that a maximum stretching of approximately 1.25 is desirable. However, the allowable 
maximum stretching may vary from case to case. In general, the best approach for a given 
case is to repeat such a grid resolution study. It should be noted that although a criterion 
for =0.1 is also employed, for the results presented in this paper, Re ce u >w = 1 proves to 
be the more strict control for these conditions. 

With the flowfield initialized to the 96-cell Cebeci-Smith results, the solution is recalcu- 
lated using the Baldwin-Lomax model. Figure 6 contrasts both numerical solutions with the 
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Note that h x is the distance from the body to the center of cell 1, so that m = 2h 1 . 

Both equations 30 and 36 provide values for An 1} and so the smallest of these two values 
is used to define the thickness of any cell adjacent to the body. With the proper stretching 
determined, dimensionality is returned to the distribution by 

n s _ 

nk = —n k (37) 

where n s is the distance from the body to the bow shock for this stack of cells (lagged from 
the previous grid). The fraction of the distance from the body to the outer grid boundary 
which lies between the bow shock and body is given by F s . A value of F s — 0.8 is used here 
so that with a converged grid, the distance across the shock layer is eighty percent of the 
total distance from the body to the outer (freestream) boundary. 


Implementation 

The turbulence models used here require a fairly developed shock layer for implementa- 
tion. Therefore, for all the cases presented herein, the LAURA code is run for laminar flow 
until the residual drops several orders of magnitude. At this point (before the laminar solu- 
tion is fully converged), the turbulence is “turned on” and the iteration process is continued. 
Not surprisingly, the switch from laminar to turbulent flow results initially in a jump in the 
residual. After a period of adjustment to the new governing equations, the residual again 
begins to steadily decrease. The rate of this decrease is less than for laminar flow, since 
the higher viscosity for turbulent flow serves as damping for the solution. Figure 4 gives a 
typical residual history. In addition to the rise with the switch from laminar to turbulent 
flow, “jumps” in the residual occur when the grid is adapted. 

The LAURA package has several features which enable it to handle very large jobs 
and utilize multiple processors when running on CRAY supercomputers. For instance, the 
computational domain can be divided into several ” blocks” with information communicated 
across their boundaries. The work within a given block (or for the entire domain) may be 
allocated to several processors. Although the results are not included here, the turbulence 
programming has been successfully tested in conjunction with both of these features. 

Solving the governing equations for turbulent flow increases the required computational 
time by approximately 1.5 percent per iteration per grid cell as compared with laminar 
calculations. However, since the turbulent terms increase the damping, more iterations are 
required to drive a turbulent solution to the same level of convergence as a laminar solution. 
Storage requirements for the LAURA code are also increased slightly since values for the 
eddy viscosity and normal distance from the body must be saved for each cell in the grid. In 
addition, the distance along the body surface to each wall-bounded cell is required to define 
the transition region. 


Results and Discussion 

In an effort to verify proper implementation of the turbulence models within LAURA, 
several test cases were computed. The results are given below in order of increasing complex- 
ity. The first case is an axisymmetric, perfect-gas flowfield calculation which is compared 


near-wall region. As a result, differences in the equilibrium models are not as significant as 
they were for the laminar results. Since the same eddy viscosity model is used for all of these 
techniques, and this is the dominant viscosity in the boundary layer, the turbulent results 
from these two algorithms should be in good agreement. This is in fact the case, as shown 
in figure 7, where the turbulent results are seen to be within approximately five percent of 
one another. 

The LAURA results presented for this case are generated using 128 cells in the normal 
direction. The maximum stretching is approximately 1.1 and occurs at the end of the body. 
Although not shown, this solution agrees well (within five percent) with a solution using 96 
cells in the normal direction, whose maximum stretching is nearly 1.2. This close agreement 
betw'een the results for two grids indicates that the solution is grid-converged. 

Three-Dimensional Flow 

The third and final case is perfect-gas flow over a generic TAV configuration, which is 
depicted in figure 8. This blended wing-body (BWB) configuration is a three-dimensional 
vehicle which features nose bluntness, a flattened windward surface, leading edge chines, 
and a double compression ramp system. It is one of several models which were tested in a 
series of studies 37 designed to explore hypersonic vehicle technology and provide data for 
computational fluid dynamic (CFD) code validation. The BWB heat-transfer and pressure 
model was tested experimentally in the Calspan shock tunnel 38 for a range of Mach and 
Reynolds numbers, at angles of attack from 0- to 10 -deg. Heat-transfer measurements were 
made at over 100 surface locations and yielded distributions along the upper and lower sym- 
metry planes, along off-centerline rays on the lower surface, and in the crossplane direction 
at several axial stations. The extensive database generated during these tests has been used 
for comparison 38_ 40 with several state-of-the-art computer codes. Comparisons are made 
with LAURA results for laminar flows over the BWB in reference 10. In the present study, 
calculations are performed for a Mach 19.6 perfect-gas case at 0 -deg angle of attack. The 
freestream conditions are p^ = 0.0066 kg/m 3 and = 240.9 K, with the wall temperature 
specified to be T w = 1256 K. 

In reference 10, the surface of the vehicle was modeled using Cheatwood and DeJar- 
nette’s ASTUD (Advanced Surface-fitting Technique with User-friendly Development) inter- 
active computer code 41 > 42 . Since the calculations were to be terminated at the end of the 
compression ramp system ( z/L — 0.778), only that portion of the geometry was modeled. 
Thus, the wings and vertical tail shown in figure 8 were not included in the surface model 
or subsequent computations. The 3-D volume grid is simply a collocation of orthogonal 2-D 
grids which were generated around cross sections of the BWB at 48 specified axial loca- 
tions. In the circumferential direction, 64 cells were used. The resultant surface distribution 
is shown in figure 9. Note that the axial spacing clusters stations around the ramp junc- 
tures {z/L = 0.435 and z/L — 0.560). For each cross-sectional plane, the circumferential 
distribution of the cells is tied to the local body curvature. 

Rather than beginning the computations for the present case from freestream values, a 
converged flowfield solution (and grid) from reference 10 is used as the initialization for the 
present case. As with the previous cases, the flowfield is solved for laminar flow initially, 
before the turbulence is “turned-on”. The transition point is specified to be s/L = .08. 
Initially, 64 cells are used in the normal direction and this converged solution supplies the 
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initialization for an 80-cell solution. The heating rates from the two solutions are generally 
within five percent of one another, so another increase in the cell density in the normal 
direction is deemed unnecessary. 

In addition to the body-surface cell distribution, figure 9 shows the converged flowfield 
grid (80-cell) at two of the cross-sectional planes. For visualization purposes, only the bound- 
aries of every fourth cell in the normal direction are shown. It should be noted that, as for 
the equilibrium case, grid adjustment proved to be cumbersome. 

As with the previous cases, the normalized heating rates (using both the Cebeci-Smith 
and Baldwin-Lomax models) are plotted as a function of nondimensional axial distance 
(where here, the axial distance is measured from virtual apex of the geometry rather than 
its blunted nose). Results are presented for the upper and lower symmetry plane, as well as 
three “off-centerline” rays on the lower surface. For selected axial stations, circumferential 
heating rates are plotted as a function of spanwise distance from the body centerline. Note 
that laminar heating-rate calculations are also presented for this case. Figure 10 provides a 
planform view of the lower surface of the body which shows the locations of the cross-sections 
and off-centerline rays where heat-transfer calculations are compared with experimental val- 
ues. 

Figure 11 shows the symmetry plane heating results for the upper and lower surface. 
Along the upper surface, the flow appears to remain laminar for the length of the body, 
although the rise in the experimental data at the end would seem to indicate the reattachment 
of a separated flow. On the lower surface, the flow is turbulent and both models predict 
the proper trends in the heating. However, the numerical results are consistently twenty to 
thirty percent higher than the data. The specified transition point may be premature, but 
that still doesn t explain the behavior for the length of the body. Similar results were seen 
in reference 10 for laminar flow along the lower centerline. 

Figure 12 gives results along three off-centerline rays. Judging from the comparison 
with the data, the flow appears to remain laminar over the forebody in the regions away 
from the centerline, with a transition to turbulence occurring at the start of the first ramp. 
In contrast to the lower centerline results, the turbulent predictions along the off-centerline 
rays are in excellent agreement with the data (generally within ten percent). In reference 10 
as well, agreement along these rays was better than along the centerline itself. 

Turning to the circumferential heating distributions, figure 13 shows the heating near 
the nose where the flow is still laminar. As a result, the three curves are nearly identical 
and are generally within fifteen percent of the data. Figure 14 shows the heating rates at 
a station just after the specified transition point. The turbulent predictions are beginning 
their departure from the laminar solution. Agreement with the data is generally within ten 
percent, except for the lower symmetry plane as noted earlier. 

Continuing downstream, figure 15 shows the heating-rate calculations at the next sta- 
tion where experimental data are available. It appears that the specified transition point is 
premature since the data more closely match the laminar solution at this station. In fact, 
excluding the upper symmetry plane, agreement between the laminar solution and the mea- 
surements is typically within ten percent. Values for the heating rates in figure 16 indicate 
the transition to turbulence has begun. Here, the lower centerline appears to be almost fully 
turbulent, although the majority of the flow at this station remains laminar. 

Recall that this case is at 0-deg angle-of- attack. As a result, because of the shape of 
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the body, the upper surface heating has been higher than the lower surface heating at the 
stations considered thus far. Figure 17 shows results at a station within the “crossover” 
region, where the lower surface heating begins to exceed the upper surface heating. Still, 
beyond a small region adjacent to the lower centerline of the body, the flow appears to 
remain laminar. The turbulent prediction is within 25 percent of the measured value, while 
differences between the laminar results and the other data range from 15 percent on the 
lower surface to 40 percent on the upper. 

The results shown in figure 18 are from a station located on the first ramp. It now ap- 
pears that the entire lower surface has transitioned to turbulent flow, perhaps triggered by 
the surface discontinuity of the forebody-ramp juncture. Agreement between the predicted 
turbulent values and the experimental results is within ten percent except along the center- 
line. The upper surface appears to remain laminar, with disagreement between calculated 
and measured values ranging from 0 to 30 percent. The numerical solution does seem to 
model the general behavior quite well, and the possibility of scatter in the experimental data 
should not be dismissed. Figure 19 shows results for a station located on the second ramp. 
Again, the lower surface appears to be fully turbulent, and agreement between calculations 
and data ranges from five to fifteen percent. The upper surface is still laminar, although the 
rise in heating on the centerline is indicative of the reattachment of a separated flow. 

The final station for comparison is presented in figure 20. As before, turbulent predic- 
tions agree well with the measured values on the ramp. The flow appears to be in transition 
in the area around the chine. On the upper surface, the separation/recirculation region is 
larger than at the previous station. 

The shock-layer flowfield can be more easily visualized in figure 21, which is a flooded 
plot of the crossflow velocity for this same station. A region of freestream flow (zero crossflow) 
lies between the outer grid boundary and the captured bow shock. A region of outflow 
extends from the bow shock down into the boundary layer. The region of darkest shading 
corresponds to the region of inflow (the boundary between these two regions is denoted as 
the dividing streamline). For both the windward and leeward sides of the body, this region 
is indicative of circumferential growth in the boundary-layer thickness. With increasing 
boundary-layer thickness, the heating decreases because the magnitude of near-wall stresses 
is reduced. In this region of inflow on the windward surface, the calculated heating does 
decrease as seen in figures 18 through 20. The fact that the results do not match the windward 
experimental values would seem to indicate more cells in the circumferential direction are 
required to accurately predict the boundary-layer growth and subsequent heating in this 
region. 

The circumferential boundary-layer growth is more dramatic on the leeward side (see 
figure 21). Separation occurs as indicated by the region very near the body where the flow 
changes direction again. This separation point corresponds to the point of minimum heating 
in figure 20. In the reattachment region, the heating rate increases from its minimum. From 
figure 20, it would appear that the flow in this region is in transition to turbulence, since 
the experimental data lie between the laminar and turbulent numerical results. However, 
without further study, it is impossible to say whether this discrepancy is in fact due to 
transition. It could be caused by inadequate grid resolution, use of the thin-layer rather 
than full Navier-Stokes equations, or shortcomings of the algebraic model as employed here. 

Figure 21 offers an opportunity to elaborate on the shortcomings of the grid alignment 
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Figure 4. Typical residual history for LAURA solution. 
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Figure 5. Grid-refinement study for 8-deg sphere-cone, R nose = 6.35cm (2.5m). 
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Figure 4. Typical residual history for LAURA solution. 
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Figure 5. Grid-refinement study for 8 -deg sphere-cone, R noS e = 6.35 cm (2.5 in). 
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Figure 6. Heat-transfer results for 8 -deg sphere-cone, R nose = 6.35 cm (2.5 in) 
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Figure 7. Heat-transfer results for 5 -deg sphere-cone, R n0 ae = 3.81cm (1.5m) 






Figure 9. Flowfield and body-surface grids for blended wing-body. 



Figure 10 . Lower planform view of blended wing-body. 
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Figure 11 . Symmetry plane heat-transfer results for the blended wing-body. 
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Figure 11. Symmetry plane heat-transfer results for the blended wing-body. 
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Figure 12. Off-centerline heat-transfer results for the blended wing-body. 
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Figure 15, Crossplane ( zjL — 0.132 ) heat-transfer results for the blended wing-body. 
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Figure 16. Crossplane { zj L = 0.179 ) heat-transfer results for the blended wing-body. 
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Figure 17. Crossplane ( zj L = 0.389 ) heat-transfer results for the blended wing-body. 
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Figure 18. Crossplane ( zjL — 0.513 ) heat-transfer results for the blended wing-body. 
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Figure 19. Crossplane ( z/L = 0.633 ) heat-transfer results for the blended wing-body. 
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Figure 20. Crossplane ( z/L = 0.743 ) heat-transfer results for the blended wing-body. 
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Figure 17* Crossplane ( z/L = 0.389 ) heat-transfer results for the blended wing-body. 
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Figure 18. Crossplane ( z/L = 0.513 ) heat-transfer results for the blended wing-body. 
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